home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / gems / g_gemsi.lha / GraphicsGems / SquareRoot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-10  |  2.0 KB  |  88 lines
  1. /*
  2.  * A High Speed, Low Precision Square Root
  3.  * by Paul Lalonde and Robert Dawson
  4.  * from "Graphics Gems", Academic Press, 1990
  5.  */
  6.  
  7.  
  8. /*
  9.  * SPARC implementation of a fast square root by table 
  10.  * lookup.
  11.  * SPARC floating point format is as follows:
  12.  *
  13.  * BIT 31     30     23     22     0
  14.  *     sign    exponent    mantissa
  15.  */
  16.  
  17. #include <math.h>
  18.  
  19. static short sqrttab[0x100];    /* declare table of square roots */
  20.  
  21. void
  22. build_table()
  23. {
  24.     unsigned short i;
  25.     float f;
  26.     unsigned int *fi = &f;    /* to access the bits of a float in  */
  27.                 /* C quickly we must misuse pointers */
  28.  
  29.     for(i=0; i<= 0x7f; i++) {
  30.         *fi = 0;
  31.  
  32.         /*
  33.          * Build a float with the bit pattern i as mantissa
  34.          * and an exponent of 0, stored as 127
  35.            */
  36.  
  37.         *fi = (i << 16) | (127 << 23);
  38.         f = sqrt(f);
  39.  
  40.         /*
  41.          * Take the square root then strip the first 7 bits of
  42.          * the mantissa into the table
  43.          */
  44.  
  45.         sqrttab[i] = (*fi & 0x7fffff) >> 16;
  46.  
  47.         /*
  48.          * Repeat the process, this time with an exponent of
  49.          * 1, stored as 128
  50.          */
  51.  
  52.         *fi = 0;
  53.         *fi = (i << 16) | (128 << 23);
  54.         f = sqrt(f);
  55.         sqrttab[i+0x80] = (*fi & 0x7fffff) >> 16;
  56.     }
  57. }
  58.  
  59. /*
  60.  * fsqrt - fast square root by table lookup
  61.  */
  62. float
  63. fsqrt(float n)
  64. {
  65.     unsigned int *num = &n;    /* to access the bits of a float in C
  66.                  * we must misuse pointers */
  67.                             
  68.     short e;        /* the exponent */
  69.     if (n == 0) return (0);    /* check for square root of 0 */
  70.     e = (*num >> 23) - 127;    /* get the exponent - on a SPARC the */
  71.                 /* exponent is stored with 127 added */
  72.     *num & = 0x7fffff;    /* leave only the mantissa */
  73.     if (e & 0x01) *num | = 0x800000;
  74.                 /* the exponent is odd so we have to */
  75.                 /* look it up in the second half of  */
  76.                 /* the lookup table, so we set the   */
  77.                 /* high bit                   */
  78.     e >>= 1;        /* divide the exponent by two */
  79.                 /* note that in C the shift */
  80.                 /* operators are sign preserving */
  81.                 /* for signed operands */
  82.     /* Do the table lookup, based on the quaternary mantissa,
  83.      * then reconstruct the result back into a float
  84.      */
  85.     *num = ((sqrttab[*num >> 16]) << 16) | ((e + 127) << 23);
  86.     return(n);
  87. }
  88.